Anomalous diffusion of proteins in sheared lipid membranes 
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We use coarse grained molecular dynamics simulations to investigate diffusion properties of 
sheared lipid membranes with embedded transmembrane proteins. In membranes without proteins, 
we find normal in-plane diffusion of lipids in all flow conditions. Protein embedded membranes 
behave quite differently: by imposing a simple shear flow and sliding the monolayers of the mem- 
brane over each other, the motion of protein clusters becomes strongly superdiffusive in the shear 
direction. In such a circumstance, subdiffusion regime is predominant perpendicular to the flow. We 
show that superdiffusion is a result of accelerated chaotic motions of protein-lipid complexes within 
the membrane voids, which are generated by hydrophobic mismatch or the transport of lipids by 
proteins. 

PACS numbers: 87.16.dj, 87.15. Vv, 87.14.cp, 82.20.Wt 
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Lipid bilayers constitute the main body of the cell 
membrane while being found in different organelles in- 
side the cell. The cell membrane hosts collections of pro- 
teins and lipid rafts and it is crowded with a variety of 
biomolccules. In such non-homogenous and diverse envi- 
ronments, the diffusion of protein molecules in lipid bilay- 
ers plays a vital role in different biological processes like 
cell signaling. The diffusion of lipids and proteins is not 
a distinct phenomenon and depends on the environment 
and neighboring molecules [l[ , and even changes from cell 
to cell |2j. Transmembrane proteins diffuse as dynamic 
complexes with lipids [H, 0] , and their interactions with 
lipid molecules mediates traffic in cell membranes. Ex- 
periments show that the hydrophobic mismatch between 
proteins and lipids controls the diffusion coefficient of 
molecules inside a bilayer [f| H| . 

Anomalous sub- and super-diffusion processes are more 
efficient scenarios for finding a nearby target than nor- 
mal diffusion 0, Q , and they enhance the formation of 
protein complexes and signal propagation. According 
to experiments, the mean square displacement (MSD) 
of membrane channel proteins of human kidney cell ex- 
hibits subdiffusio n 3 . The addition of cholestcrols to 
lipid membranes [lCj . and the augmented area coverage 
of membrane proteins [ll[ also lead to subdiffusion of 
lipids and proteins. Superdiffusivity has been observed 
in several physical systems [TU, [l3| . A recent study by 
Kohlcr et al. [l4[ shows that in a gel composed of actin 
filaments, fascin molecules and myosin-II filaments, the 
diffusion of small actin and fascin clusters are superdif- 
fusive because of the work done by molecular motors. 

In many conditions membranes are under shear. When 
a red blood cell (RBC) migrates through vessels smaller 
in diameter than itself, the RBC membrane is under 
shear. The blood flow exerts tangential shear stresses on 
vascular endothclia, and initiates cellular processes like 
activating G protein-coupled receptors. These receptors 
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are able to sense the fluid shear stress as an increase in the 
lateral membrane tension, and subsequently go through 
conformational changes [lfj. The temporal and spatial 
changes in the membrane fluidity, in response to shear 
flow, have been observed experimentally [3, E3] ■ 

In this study we are interested in the diffusivity of lipid 
and protein molecules in flat membranes under shear 
flow, and attempt to answer three fundamental questions 
using molecular dynamics (MD) simulations: (i) Do lipid 
and protein molecules have different diffusion properties 
parallel and perpendicular to flow direction? (ii) How 
does a simple shear flow influence the random motions of 
transmembrane proteins? 

Model and methods. We simulate lipid membranes uti- 
lizing a flexible lipid model and triple-strand rigid 
proteins flS fl [F ig. Ufa)]. We use the model of Khosh- 
nood et al. [19| and perform MD simulations of an NVT 
ensemble. Although different coarse grained models have 
been developed over years [2(|, the model adopted here 
has the ability to mimic the physical properties of lipid 
membranes. We express the position and velocity vectors 
of particles in the Cartesian [x, y,z) coordinate system 
whose origin is located at the center of our cubic simu- 
lation box. The x and y axes lie in the membrane plane 
and the z axis is perpendicular to that. 

MD scales of length, time, mass and energy are 
(j = 1/3 nm, t = 1.4 ps, N avo m = 36 g/mol and 
N avo e = 2 kJ/mol, respectively. N avo is the Avogadro's 
number. In all simulations, the dimensions of the box 
along the coordinates axes, L x , L y and L z , are set to 
L x = L y = L z = 28.71cr. The total number of particles 
equals N = 15625, which gives a fixed number density 
p = 2/(3er 3 ). Lees-Edwards boundary condition is em- 
ployed to generate simple shear flow [2l| with the shear 
rate 7 = 0.03r _1 . Other boundary conditions are peri- 
odic. The temperature is set to 324 K so that the system 
is safely above the gel to liquid phase transition tempera- 
ture of different phosphatidylcholine lipid bilayers. Phys- 
ical quantities are measured using a run-time of ~ 5000r. 
We compute the tension of our model membranes from 



{Pxx + Pyy) /2]L Z where P aa (a = x,y,z) 
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FIG. 1. (a) The models of lipid and protein molecules. Red 
and yellow spheres are hydrophilic and hydrophobic parti- 
cles, respectively, (b) Velocity profile of a sample solvent- 
membrane system under simple shear flow with 7 = 0.03r _1 . 

are the components of the pressure tensor. 

In this study we apply MSD to determine the diffusion 
properties of randomly moving particles. The diffusion 
coefficient is thus calculated using the Einstein expression 

1 N 

D aa = Hm — [©«(*) - MO)] 2 }, (1) 

i=l 

where a £ {x, y, z} and qi a is the displacement due to 
the random motion of the ith particle in the a-direction. 
The summation in Eq. ([1]) is taken over the particles 
of the same type. From here on, we will drop the sum- 
mation sign for brevity. The operator (• • ■ } denotes the 
canonical average. Eq. (JlJ describes the regular Brown- 
ian motion when the MSD is linearly proportional to t. 
For anomalous diffusion one has 

{[q a (t)-q a (0)f)=2D a a J a , (2) 

where a is the diffusion exponent and D^ a is the frac- 
tional diffusion coefficient. The regimes with < a < 1 
and a > 1 are subdiffusive and superdiffusive, respec- 
tively. To obtain smooth MSD curves, we evolve sys- 
tems of 80 different initial conditions and report their 
ensemble-averaged diffusion coefficients and MSDs. The 
diffusion of lipids is investigated by tracing the motion of 
their head groups. Proteins are traced using their center 
of masses. 

In equilibrium models the flux of particles is associ- 
ated with their random motions, and any displacement 
is due to thermal fluctuations. In sheared membranes, 
however, there is a combination of streaming and diffu- 
sive fluxes. We thus need to distinguish and eliminate 
the streaming flux when calculating the MSD. Let us de- 
fine the actual velocity components of the jth particle as 
Vj a = (v a ) +Vja where (v a ) is the average streaming ve- 
locity, and Vj a is the peculiar velocity whose time integral 
gives the displacements in Eqs. ([T]) and ([2]). In equilib- 
rium models, the average velocities (v a ) vanish and we 
obtain qj a = J Vj a dt = J Vj a dt. With external shearing, 
the flow is always imposed in the x-direction. Therefore, 



Vj y = Vjy and Vj Z = Vj z are directly integrated to find the 
corresponding displacements. When the simulation box 
is uniformly filled with one type of particles (let us says 
solvent particles), one readily finds Vj X = Vj X — zj. In 
the presence of a lipid bilayer, the vertical velocity pro- 
file in the z-direction is no longer linear [see Fig. [IJb)]. 
Therefore, to obtain the MSD of lipids, we define (v x ) as 
the average velocity of the layer where the head groups of 
phospholipds reside, and obtain qj X = J (vj x — (v x )) dt. 

Diffusion of lipids. The lateral diffusion of lipids in 
equilibrium conditions is enhanced as the membrane ten- 
sion increases [22|, H3|- By turning on the shear flow, 
lipid molecules undergo an initial ballistic motion that 
transforms into an interval of subdiffusion with a = 0.7 
[Fig. Ufa)]. The transient anomalous state has been ob- 
served in atomistic simulations [24| as well. After the 
transient anomalous diffusion and over longer time scales 
a normal diffusion with a = 1 is observed [Fig. HJa)]. 
It is noted that we have found similar MSD profiles for 
lipid molecules in equilibrium and sheared systems, and 
in both cases lipids ultimately develop normal diffusion. 
Although Kneller et al. [25[ reported a permanent sub- 
diffusive behavior by the atomistic simulation of lipid 
membranes in equilibrium, experiments support a final 
regular diffusion regime, as we do, even in the presence 
of obstacles (2|| . We conclude that the diffusion regime 
of lipid molecules is invariant with and without external 
shearing. 

The diffusion coefficients obtained from the normal dif- 
fusion region of MSD plots, are larger for smaller shear 
rates. For example, for a membrane of Ni = 600 lipid 
molecules, we find C = (1-4846 ± 0.2624)cr 2 /e, D xx = 
0.001608o- 2 /t and D yy = 0.001587ct 2 /t. For the same 
system under a shear flow of 7 = 0.03r _1 , the mem- 
brane tension drops to f = (0.8163 ± 0.2727)er 2 /e and 
diffusion coefficients reduce to D xx = 0.001513ct 2 /t and 
Dyy = 0.001580(T 2 /r. The reason is that the membrane 
thickness increases for higher shear rates and the tension 
decreases without any change in the area per lipid [l9l ]. 
Consequently, the fluidity of the membrane decreases and 
slows down the diffusion process. After applying the 
shear force, we find that D xx drops for about 6 percent 
while Dyy remains almost constant [this is indistinguish- 
able in Fig. Ufa)]. We speculate that the alignment of 
lipid chains with the flow breaks the isotropy and yields 
D xx 7^ Dyy. In the z-dircction, perpendicular to the 
membrane plane, our MSD plots always show a confined 
motion as is expected. 

Diffusion of proteins. We add rod-like proteins to the 
membrane, and simulate models with different protein 
concentrations that vary significantly from cell to cell. 
Since proteins increase the membrane tension as they 
perturb the distribution of lipids [l9j . we increase the 
number of lipids (proportional to proteins) to keep the 
membrane tensionless. In equilibrium and for a mem- 
brane with a single embedded protein, we can mea- 
sure the diffusion coefficients D aa since the MSD of 
protein shows an ultimate regular diffusion. We find 
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FIG. 2. (a) MSD of lipid molecules. The membrane is sheared 
and has 600 lipids, (b) MSD of protein molecules in a protein- 
embedded sheared model with 2 and 4 proteins, which corre- 
spond to 640 and 660 lipids, respectively. 
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yy — 0.0254cr 2 /r, which are 
equivalent to D xx k, D yy pa 2 x 10~ 9 m 2 /s. These val- 
ues are larger than experimental values, by two orders 
of magnitude. The obvious reason is the effect of coarse 
graining that has reduced the interdigitation and fric- 
tion between molecules, and allows for faster movements 
of particles. By putting the system under simple shear 
flow, proteins undergo Brownian motion when only two 
proteins are used [Fig. [2jb)]. One could anticipate this 
result, for single proteins cannot remarkably perturb the 
distribution of lipids, and change the diffusion proper- 
ties of the membrane. With 4 proteins, however, we ob- 
serve that they form two double-protein clusters (due to 
the depletion force 1 ), and exhibit a strong superdiffusive 
motion parallel to the flow. Fig. &b) shows how af- 
ter t ~ lOOr the normal diffusion regime transforms to 
strong superdiffusion with a = 1.7 in the x-direction. 
Interestingly, this exponent is the same as the superdif- 
fusion exponent found by Kohler et al. 14[ for active 
diffusion of protein clusters by molecular motors. Be- 
cause of crowding effect and increase in the concentra- 
tion of proteins our results show a subdiffusive be- 
havior along the y axis with a = 0.7. Weigel et al. [i| 
observed a = 0.8 ± 0.1 in experiments with channel pro- 
teins of human kidney cell, and Javanainen et al. (27| 
found a ~ 0.75 ±0.15 by molecular simulations of aggre- 
gating NaK channel proteins. 

Let us define the local concentration of the head par- 
ticles of lipid and protein molecules at the position r 
and time t as / = J2?=i i H ( S i) ~ H ( S i ~ A )l where 
Si(t) = \ri(t) — r\, and ri(t) is the position vector of 
the ith head particle. Ah denotes the total number of 
head particles in the monolayer, is the Heavisidc 

step function, and 2A is the typical size of the cross sec- 
tion of a protein cluster (or a protein-lipid complex). 
Our numerical experiments show A = 4er is the best 
choice. We examine the trajectories of protein molecules 
and the spatial variation of the normalized distribution 
f(r,&,t) = (/ - / mi n)/(/max - /min) to explain the 
physics of observed superdiffusion. Here / m ; n and / max 
are the minimum and maximum values of / at a given 
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FIG. 3. The local concentration f(r,A,t) of head groups in 
the upper (left panel) and lower (right panel) leaflets of a 
sheared membrane. The head particles of proteins are shown 
by filled circles. Drawn circles have radii of A, and their 
centers lie at the centroid of the head particles of proteins. 
The flow is in the x + and x~ directions for upper and lower 
leaflets, respectively [cf. Fig. [TJb)]. The model has 660 lipids 
and two double-protein clusters. 



time t. Fig. [3] demonstrates contour plots of / for the 
upper and lower monolayers at a randomly selected time. 

Fig. @|a) demonstrates the trajectories of a single pro- 
tein molecule and two double-protein clusters, in equi- 
librium and sheared systems, respectively. The equilib- 
rium trajectory corresponds to regular diffusion because 
it covers a definite area. In the sheared system trajec- 
tories show local isotropic wanderings followed by small- 
step jumps mainly in the flow direction. These successive 
jumps can be interpreted by inspecting the contour plots 
of /(r, A, t) over a long duration of time. The hydropho- 
bic mismatch between protein clusters and the membrane 
disturbs the bilayer thickness and the arrangement of 
nearby lipids. Moreover, proteins are able to transport 
their neighboring lipids with them and behave as dy- 
namic complexes 0, These two effects collaborate 
to create transient voids whose distribution can be de- 
scribed by g = 1 — f (light shades in Fig. [3]). When the bi- 
layer is sheared, protein-lipid groups are pushed into the 
voids created by themselves or other groups/complexes 
and experience accelerated, and therefore, superdiffusive 
movements. It should be noted that during our simula- 
tions, the center of mass of the membrane and embedded 
proteins remains almost at the center of the coordinate 
system. 

We have computed the probability distribution func- 
tion (PDF) of protein displacements and plotted it in 
Fig. HJb). A Gaussian function has been fitted to the 
data by setting its maximum to the maximum of PDF, 
and its variance is found using the full width at half mini- 
mum of PDF. The PDF exhibits a deviation from normal 
distribution and it has tails. We have also applied the 
Kolmogorov-Smirnov test [28| to confirm that the PDF 
is not a normal Gaussian. This is a clear indication of 
anomalous diffusion. 

The ends of proteins are pulled in opposite directions 
by the two sheared solvent columns. An important ques- 
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FIG. 4. (a) Trajectories of a single protein in equilibrium 
(red) and two double-protein clusters (cyan and green) under 
shear flow, (b) PDF for the displacements of a sample two- 
protein cluster. Solid line shows the best Gaussian function 
fitted to the data, (c) Autocorrelation function -A(ii ag ). The 
correlation time t c ~ 23r is defined at the point where yl(ti ag ) 
abruptly drops below 10% of its maximum. We have taken 
1000 successive samples in the time domain, with increments 
of It, to compute A. (d) Cross-correlation function C(£i ag ). 
The integrals in C have been taken using a grid of 29 x 29 in 
the xj/-plane and 1000 successive points, with steps of lr, in 
the time domain. 



tion is why do protein clusters prefer to jump into the 
voids when the membrane is sheared? We have a simple 
explanation for this behavior: Two ends of proteins at- 
tract lipids from two different layers of the membrane. 
Thus, the symmetry in the distribution of the upper 
and lower protein-bound lipids is likely to break. More- 
over, the shear force is exerted by the solvent on both 
the lipid and protein heads. The mentioned symme- 
try breaking thus leads to different force components 
at the upper and lower ends of protein-lipid complexes, 
and they will jump into a void in the direction specified 
by the broken symmetry. To quantify this process, we 
take a sample two-protein cluster and define r u and ri 
as the center of masses of its protein heads in the up- 
per and lower leaflets, respectively. We then compute 
= f( r u, A,tj) - f(rt,A,ti), which is proportional 
to the net shear force exerted on the cluster at the time 



step U: the local effective area in contact with a solvent 
column is determined by the number of head particles, 
and the shear force is calculated by multiplying the ef- 
fective contact area by the shear rate and viscosity. 
will be zero if the concentrations of the head particles of 
lipids are identical around the two heads of the cluster. 
Defining £ as the average of the autocorrelation 

function 

l "" )= E, [««.)-«] " ' " 

plotted in Fig. @Ic) carries interesting information about 
the shear force experienced by the cluster: the correla- 
tion time t c « 23t shows a sustained accelerated mo- 
tion of the cluster over to < t < to + t c , indepen- 
dent of the initial time to. Moreover, the oscillatory 
decaying profile of A(ti ag ) shows a random sy mmetry 
breaking in the sense of deterministic chaos |29|, §5.3.4]. 
The existence of a correlation time can also be veri- 
fied by studying the cross-correlation function C(t\ ag ) = 
J dt J drf(r, A, t+t\ ag ) -g(r, A, t) over one of the leaflets. 
Fig. 0Jd) shows that C(t\ ag ) of the upper leaflet steeply 
rises from t\ a g = and peaks at ti ag ps 21t, which is the 
earliest time span that protein-lipid complexes need to 
occupy their nearest voids. This is quite consistent with 
the acceleration time scale of protein clusters predicted 
by the autocorrelation function A{t\ ag ). For ti ag > 21r 
the cross-correlation function remains almost flat because 
the size of a void is always bigger than the distance that 
a protein cluster travels during t ~ 0(t c ). 

Conclusions. Cell responses to stimuli arc fast due 
to enhanced mobility of protein receptors. Consider our 
findings, supcrdiffusion of proteins under shear flow can 
play a dominant role in the process of signaling in en- 
dothelium cells, RBCs, liposomes used for targeted drug 
delivery and other sheared membranes. We note that 
since the length and shape of proteins and their ability 
in attracting neighboring lipids controls the sizes of voids 
— there was no other mechanism in our models for void 
generation — supcrdiffusion properties reported in this 
study are not universal and they highly correlate with 
the properties of embedded proteins. Simulations using 
detailed structures of lipids and proteins are needed to 
better assess the superdiffusive behavior in realistic cell 
membranes. Experimental exploration of our results can 
be done using single particle tracking [30j which can give 
the MSD of proteins directly. 
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